home *** CD-ROM | disk | FTP | other *** search
/ Libris Britannia 4 / science library(b).zip / science library(b) / MATH / POLYSIM.ZIP / POLYSIM.C next >
Text File  |  1989-05-21  |  9KB  |  275 lines
  1. /* Program POLYSIM.C
  2.  
  3.    (c)   D.J. Murphy  1989                  Borland TurboC v1.5
  4.  
  5.    Room 213
  6.    Chemistry Dept.
  7.    University of Edinburgh
  8.    King's Buildings
  9.    West Mains Road
  10.    Edinburgh
  11.    Scotland
  12.    EH9 3JJ
  13.  
  14.    JANET: djm@uk.ac.etive         ARPA:djm%ed.etive@nsfnet-relay.ac.uk
  15.           Murff@uk.ac.edinburgh        Murff%ed@nsfnet-relay.ac.uk
  16.           trinity@uk.ac.ed.cs.tardis   trinity%ed.cs.tardis@nsfnet-relay.ac.uk
  17.  
  18.   Description:
  19.   Uses a modified simplex method to optimize a best-fit polynomial
  20.   function to a set of input x,y data pairs. The polynomial is of the form:
  21.  
  22.   y = ax^5 + bx^4 + cx^3 + dx^2 + ex + f + g/x + h/x^2
  23.  
  24.   and POLYSIM optimizes the coefficients a through h to minimize the mean
  25.   deviation of y(calculated) from y(given).
  26.   Note that, since the optimization is based on MEAN deviations within the
  27.   range of the given dataset, it is not recommended that this be used for
  28.   extrapolations greatly outwith the input range.
  29.  
  30.   Reference :
  31.   S.N. Deming & S.L. Morgan Anal. Chem. 45, 278A(1973)
  32.  
  33. */
  34.  
  35. /************ HEADERS **************/
  36.  
  37. #include <stdio.h>
  38. #include <process.h>
  39. #include <math.h>
  40. #include <stdlib.h>
  41.  
  42. /************ DEFINES **************/
  43.  
  44.  
  45. #define MAXDAT   100         /* Max # of data points */
  46. #define BOUNDARY 0.01        /* Boundary condition = lowest mean dev
  47.                                  BOUNDARY * y-data range */
  48.  
  49. /***********************************/
  50.  
  51. main(argc, argv)
  52.  
  53. int   argc;
  54. char  *argv[];
  55.  
  56. {
  57.  double         xinput[MAXDAT], yinput[MAXDAT], meandev[9], md2[9];
  58.  double         ycalc[MAXDAT][9], cf1[8], yrange;
  59.  double         coeff[9][8];
  60.  double         adump, bdump, cdump[MAXDAT];
  61.  int            count = 0;
  62.  int            sortdump[9];
  63.  int            acount, bcount, ccount, flag;
  64.  unsigned long  dcount = 0;
  65.  char           xin[20], yin[20];
  66.  FILE           *infile, *temp;
  67.  
  68.  /* Initialize starting coefficients */
  69.  
  70.  for (acount = 0; acount < 8; acount++) {
  71.      for (bcount = 0; bcount < 9; bcount++)
  72.             coeff[bcount][acount] = 0;
  73.  }
  74.  for (acount = 0; acount < 9; acount++) {
  75.      for (bcount = 0; bcount < acount; bcount++)
  76.      coeff[acount][bcount] = 1;
  77.  }
  78.  
  79.  /* Check command line parameters and get input file */
  80.  
  81.  if (argc != 3) {
  82.      fprintf(stdout, "Usage: polysim inputfile.ext outputfile.ext");
  83.      exit(1);
  84.  }
  85.  
  86.  if ((infile = fopen(argv[1], "r")) == NULL) {
  87.     fprintf(stdout, "Cannot open data file %s\n", argv[1]);
  88.     exit(1);
  89.  }
  90.  while (fscanf(infile, "%s %s", xin, yin) != EOF) {
  91.        xinput[count] = atof(xin);
  92.        yinput[count] = atof(yin);
  93.  
  94.  /* Vet data to get rid of values likely to lead to overflow or
  95.     divide by zero errors */
  96.  
  97.        if (fabs(xinput[count]) > 1e-150)
  98.           count++;
  99.  }
  100.  fclose(infile);
  101.  
  102.  /* Find range of y-data for termination check later on */
  103.  adump = 1e300;        /* Will become lowest y */
  104.  bdump = -1e300;       /* Will become highest y */
  105.  for (ccount = 0; ccount < count; ccount++) {
  106.      if (yinput[ccount] < adump)
  107.     adump = yinput[ccount];
  108.      if (yinput[ccount] > bdump)
  109.     bdump = yinput[ccount];
  110.  }
  111.  yrange = bdump - adump;
  112.  
  113.  /* Set up simplex */
  114.  
  115.  flag = 1;
  116.  
  117.        /* Flush ycalc */
  118.  
  119.        for (acount = 0; acount < count; acount++) {
  120.        for (bcount = 0; bcount < 9; bcount++)
  121.            ycalc[acount][bcount] = 0;
  122.        }
  123.  
  124.        /* Evaluate first ycalc */
  125.  
  126.        for (ccount = 0; ccount < count; ccount++) {        /* Data point */
  127.            for (acount = 0; acount < 9; acount++) {        /* Vertex 0 - 8 */
  128.                for (bcount = 7; bcount >= 0; bcount--)      /* Coefficient 0 - 7 */
  129.                    ycalc[ccount][acount] += (coeff[acount][bcount] * pow(xinput[ccount], (bcount - 2)));
  130.            }
  131.        }
  132.        for (acount = 0; acount < 9; acount++) {         /* Evaluate mean deviations */
  133.        meandev[acount] = 0;                            /* flush meandev */
  134.            for (bcount = 0; bcount < count; bcount++)      /* for each vertex */
  135.                meandev[acount] += (ycalc[bcount][acount] - yinput[bcount]);
  136.  
  137.            meandev[acount] = fabs((meandev[acount] / count));
  138.  
  139.        }
  140.        /* Start running simplex */
  141.  
  142. while (flag) {
  143.  
  144.        ++dcount;
  145.  
  146.        for (acount = 0; acount < 9; acount++)
  147.            md2[acount] = meandev[acount];
  148.  
  149.        /* Sort results by meandev (lowest to highest) into index array */
  150.  
  151.        for (acount = 0; acount < 9; ++acount) {
  152.            adump = 1e300;
  153.            for (bcount = 0; bcount < 9; ++bcount) {
  154.                if (md2[bcount] < adump) {
  155.               ccount = bcount;
  156.                      adump = md2[bcount];
  157.            }
  158.          }
  159.            md2[ccount] = 2e300;
  160.            sortdump[acount] = ccount;
  161.        }
  162.  
  163.        /* Flush md2 */
  164.  
  165.        for (acount = 0; acount < 9; acount++)
  166.        md2[acount] = 0;
  167.  
  168.        /* Find centroid of first 8 points and flush cf1 */
  169.        for (acount = 0; acount < 8; acount++) {
  170.        for (bcount = 0; bcount < 8; bcount++)
  171.            md2[acount] += coeff[(sortdump[bcount])][acount];
  172.  
  173.            md2[acount] /= 8;       /* md2 now contains coordinates of centroid */
  174.        cf1[acount] = 0;
  175.        }
  176.  
  177.        /* Then reflect last point through it and flush cdump */
  178.  
  179.        for (acount = 0; acount < 8; acount++) {
  180.        cf1[acount] = 2 * md2[acount] - coeff[(sortdump[8])][acount];
  181.            cdump[acount] = 0;
  182.        }
  183.  
  184.        /* Then evaluate response */
  185.  
  186.        for (acount = 0; acount < count; acount++) {
  187.            bdump = 0;
  188.        for (bcount = 7; bcount >= 0; bcount--)
  189.            bdump += cf1[bcount] * pow(xinput[acount], (bcount - 2));
  190.  
  191.            cdump[acount] += (bdump - yinput[acount]);
  192.        }
  193.        for (acount = 0; acount < count; acount++)
  194.        bdump += cdump[acount];
  195.        bdump = fabs((bdump / count));  /* bdump now = response of new point */
  196.  
  197.  
  198.        /* See how new response compares with what we already have and act
  199.       accordingly:
  200.       if bdump < meandev[(sortdump[0])] then double point translation
  201.       if bdump > meandev[(sortdump[8])] then quarter point translation
  202.       if bdump > meandev[(sortdump[0])] && bdump < meandev[(sortdump[4])]
  203.           then 75% translation
  204.       otherwise keep translation.
  205.  
  206.           Variable status:
  207.       coeff .... coefficients from last cycle
  208.       cf1 ...... coefficients from simple reflection of worst point
  209.       meandev .. mean deviations from last cycle
  210.       dbump .... mean deviation of cf1
  211.       sortdump . index to meandev & coeff of quality of solutions sorted
  212.              lowest to highest in meandev
  213.       md2 ...... coordinates of centroid of best 8 functions
  214.        */
  215.  
  216.        if (bdump < meandev[(sortdump[0])]) {
  217.       for (acount = 0; acount < 8; acount++)
  218.           coeff[(sortdump[8])][acount] = 4 * md2[acount] - 2 * coeff[(sortdump[8])][acount];
  219.        } else if (bdump > meandev[(sortdump[8])]) {
  220.          for (acount = 0; acount < 8; acount++)
  221.              coeff[(sortdump[8])][acount] = ((md2[acount]/2) - (coeff[(sortdump[8])][acount]/4));
  222.        } else if ((bdump > meandev[(sortdump[0])]) && (bdump < meandev[(sortdump[4])])) {
  223.          for (acount = 0; acount < 8; acount++)
  224.              coeff[(sortdump[8])][acount] = ((6 * md2[acount] - 3 * coeff[(sortdump[8])][acount]) / 4);
  225.        } else {
  226.               for (acount = 0; acount < 8; acount++)
  227.               coeff[(sortdump[8])][acount] = cf1[acount];
  228.        }
  229.  
  230.        /* Recalculate point */
  231.  
  232.        for (acount = 0; acount < count; acount++) {
  233.            bdump = 0;
  234.        for (bcount = 7; bcount >= 0; bcount--)
  235.            bdump += coeff[(sortdump[8])][bcount] * pow(xinput[acount], (bcount - 2));
  236.  
  237.            meandev[(sortdump[8])] += (bdump - yinput[acount]);
  238.        }
  239.        meandev[(sortdump[8])] = fabs((meandev[(sortdump[8])] / count));
  240.  
  241.        /* Check for termination */
  242.  
  243.        adump = fabs((meandev[(sortdump[0])]/yrange));
  244.        if (adump < BOUNDARY)
  245.       flag = 0;
  246.  
  247.        /* Loop on flag = 1 */
  248.  
  249.  }
  250.  
  251.  /* report results */
  252.  
  253.  if ((infile = fopen(argv[2], "w")) == NULL) {
  254.     fprintf(stdout, "\nCannot write to output file %s", argv[2]);
  255.     exit(1);
  256.  }
  257.  fprintf(infile, "Results of analysis of data file : %s after %u iterations\n", argv[1], dcount);
  258.  fprintf(infile, "\nCoefficients :");
  259.  for (acount = 0; acount < 8; acount++) {
  260.      if (fabs(coeff[0][acount]) < 1e-20)
  261.     coeff[0][acount] = 0;
  262.  
  263.      fprintf(infile, "\n%c %15g", (acount + 97), coeff[0][acount]);
  264.  }
  265.  fprintf(infile, "\n\nMean deviation = %15g", meandev[(sortdump[0])]);
  266.  fprintf(infile, "\n\nPoint         X-given         Y-given    Y-calculated      Difference");
  267.  for (acount = 0; acount < count; acount++) {
  268.      adump = fabs((yinput[acount] - ycalc[acount][(sortdump[0])]));
  269.      fprintf(infile, "\n%5d %15g %15g %15g %15g", (acount + 1), xinput[acount], yinput[acount], ycalc[acount][(sortdump[0])], adump);
  270.  }
  271.  
  272.  fclose(infile);
  273.  
  274. }
  275.